MODULE md1_shell_pairs

   USE md1_global_types
   USE md1_global_consts
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_get_shell_pairs,    &
             md1_free_shell_pairs

   TYPE(md1_sh_pairs), ALLOCATABLE, TARGET :: sh_pairs(:)

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_get_shell_pairs(basis,sh_pairs_ptr)

      TYPE(md1_basis), INTENT(IN) :: basis
      TYPE(md1_sh_pairs), POINTER :: sh_pairs_ptr(:)
      INTEGER(INTK) :: n_pairs

      IF (ALLOCATED(sh_pairs)) THEN
         sh_pairs_ptr => sh_pairs
      ELSE
         ! make list of non-vanishing shell pairs
         CALL md1_make_shell_pairs(basis,n_pairs)  ! first get n_pairs
         ALLOCATE(sh_pairs(n_pairs))
         CALL md1_make_shell_pairs(basis,n_pairs)  ! now store list
         sh_pairs_ptr => sh_pairs
         WRITE(IO1,*) 'Number of shells      =', basis%nshell
         WRITE(IO1,*) 'Number of shell pairs =', SIZE(sh_pairs)
      END IF

   END SUBROUTINE md1_get_shell_pairs

!-------------------------------------------------------------------------------

   SUBROUTINE md1_free_shell_pairs

      IF (ALLOCATED(sh_pairs)) DEALLOCATE(sh_pairs)

   END SUBROUTINE md1_free_shell_pairs

!-------------------------------------------------------------------------------
!
!   FUNCTION fmm_extent(ExpPI)
!
!!      USE fmm_md4_globals, ONLY: fmm_shell_pair_thresh, fmm_X0
!!      USE fmm_md4_globals, ONLY: fmm_grain_inv, fmm_extent_min
!!      USE fmm_box_utils,   ONLY: fmm_branch
!
!      IMPLICIT NONE
!
!      REAL(REALK), INTENT(IN) :: ExpPI
!
!      REAL(REALK), PARAMETER :: root3 = SQRT(3.0d0)
!      REAL(REALK), PARAMETER :: thresh = 10d0  ! 1e-10 screening
!      REAL(REALK)   :: tmp1, tmp2
!      REAL(REALK)   :: fmm_extent
!      INTEGER(INTK) :: bra_min
!
!!      REAL(REALK), PARAMETER :: erfc_inv = 3.4589  ! 1e-6
!!      REAL(REALK), PARAMETER :: erfc_inv = 4.0522  ! 1e-8
!!      REAL(REALK), PARAMETER :: erfc_inv = 4.5728  ! 1e-10
!!      REAL(REALK), PARAMETER :: erfc_inv = 5.6739  ! 1e-15
!
!!      fmm_extent = fmm_X0*ExpPI
!
!      ! Extent based on Classical overlap
!!      tmp1 = SQRT(ExpPI)*erfc_inv
!
!!      ! Artificial extent to keep all near-field exact integrals
!!      bra_min = fmm_branch(fmm_extent_min,fmm_grain_inv)
!!      tmp2 = half*root3*(bra_min+1)/fmm_grain_inv
!
!!      tmp2 = zero
!!      fmm_extent = MAX(tmp1,tmp2)
!
!      fmm_extent = SQRT(thresh*ExpPI*(LOG_10))
!
!   END FUNCTION fmm_extent
!
!-------------------------------------------------------------------------------

   SUBROUTINE md1_make_shell_pairs(basis,np)

      IMPLICIT NONE

      TYPE(md1_basis), INTENT(IN)  :: basis
      INTEGER(INTK),   INTENT(OUT) :: np

!fixme
      REAL(REALK), PARAMETER :: thresh = -2.30258d+00 * 13.0d0
      TYPE(md1_sh_pairs), POINTER :: sh
      INTEGER(INTK) :: Ish,Jsh, I,J,Jhi, IJ
      REAL(REALK)   :: A(3), B(3), M(3), P(3), PM(3), RAB(3)
      REAL(REALK)   :: ExpA, ExpB, ExpP, ExpPI, ExpAR2
      REAL(REALK)   :: R2AB, ExpKAB

      np = 0
      Ishel: DO Ish = 1, basis%nshell
         A(:) = basis%xyz(:,basis%atom(Ish))
         Jshel: DO Jsh = 1, Ish
            B(:) = basis%xyz(:,basis%atom(Jsh))

            M(:) = Half * (A(:) + B(:))
            RAB(:) = A(:) - B(:)
            R2AB = DOT_PRODUCT(RAB,RAB)
      
            IJ = 0
            DO I = basis%exp_id(Ish), basis%exp_id(Ish) + basis%nprim(Ish) - 1
               ExpA = basis%exp(I)
               ExpAR2 = ExpA * R2AB

               Jhi = basis%exp_id(Jsh) + basis%nprim(Jsh) - 1
               IF (Ish == Jsh) Jhi = I
               DO J = basis%exp_id(Jsh), Jhi
                  ExpB = basis%exp(J)
                  ExpP = ExpA + ExpB
                  ExpPI = One / ExpP
                  ExpKAB = - ExpAR2 * ExpB * ExpPI

                  IF (ExpKAB >= thresh) THEN
                     IJ = IJ +1
                     P(:) = (ExpA * A(:) + ExpB * B(:)) * ExpPI
                     PM(:) = P(:) - M(:)
                  END IF

               END DO
            END DO
      
            IF (IJ <= 0) CYCLE
            np = np + 1
            IF (.NOT. ALLOCATED(sh_pairs)) CYCLE

            IF (np > SIZE(sh_pairs)) STOP 'fixme'
            sh => sh_pairs(np)
            sh%I = Ish
            sh%J = Jsh
            sh%xyz(:) = M(:)
            sh%Iang = basis%L(Ish)
            sh%Jang = basis%L(Jsh)
            sh%P%lo = basis%car_AO(Ish)
            sh%Q%lo = basis%car_AO(Jsh)
            sh%P%hi = basis%car_AO(Ish) + basis%LM(basis%L(Ish))%car_nlm -1
            sh%Q%hi = basis%car_AO(Jsh) + basis%LM(basis%L(Jsh))%car_nlm -1
   
         END DO Jshel
      END DO Ishel

   END SUBROUTINE md1_make_shell_pairs

!-------------------------------------------------------------------------------

END MODULE md1_shell_pairs
